% Initialises Spectral Decomposition
% Christoph Gortz
%
% Calls: - get_all_parameter_values
%        - generate_model_structures
%        - perform_spectral_decomp, calls codes by Alejandro Justiniano
%
% Notes: - dynare file xxx_restults.m has to be in the directory! 
%          Call the file estimation_results.m!
%          Name has to be specified below. 
%          Name has to differ from xxx_results.m generated by .mod file
%          otherwise it will be overwritten later!
%        - specify name of dynare file to call below. This .mod file has to
%          be in the same directory.
%        - make specifications in all sections of perform_spectral_decomp!!
%        



% CHANGES REGARDING THE CORRESPONDING DYNARE FILE!!!!!!!!!
% dynarefile.mod has to be in the directory! Name has to be specified below.
%
% The dynare-model file and the _results.m file from the estimation need to
% have different names otherwise the latter one will be overwritten.
%
% Include the following lines at the top in dynarefile.mod
%      // Load calibration of the estimated parameter values
%      load ParameterValuesForDynare
%
% also delete all specified parameter values in the .mod file of parameters
% that are estimated. Replace them with e_parametername - this info is
% loaded from DecompInfo
%
% also specify something similar to "stoch_simul(irf = 0, periods = 360, noprint)"
%
% In the shocks-section, replace the numerical value for the shock
% variances by std_shockname.
%
% Comment out the "steady" command in the dynare .mod file. Steady state
% values will not be printed which makes routine faster.
%
% Variable desc has to be smaller than simulation periods specified in .mod file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clc
clear


%% INITIALISATION

calibmethod = 1;  % Determines how parameter values are chosen
                  % If calibmethod == 1, parameter values are chosen
                  % according to the posterior densities. 
                  % If calibmethod == 2, parameter 
                  % values are the posterior mode.

L = 50;            % Number of Loops per shock specification
if calibmethod == 2;
    L = 1;
end

dynareFileName = 'twosb_v12b_B_TFP48restrrho'; % Name of dynare file to be called

desc = 100;       % Descarded periods at the beginning of the time-series
lp = 5;           % Lower percentile (e.g. 10 for 10% percentile)
up = 95;          % Upper percentile (e.g. 90 for 10% percentile)

% Define number of Shocks
load estimation_twosb_v12b_B_TFP48restrrho_results   % Contains posterior densities of 
%                           estimated parameters
Postdens2 = oo_ .posterior_density.shocks_std;
numshocks = numel(fieldnames(Postdens2));
Ns = numshocks;           % Number of Shocks

save Store_permanent L Ns desc lp up calibmethod 


%% GET PARAMETER VALUES
% Assigns parameter values from posterior densities or mode
for l = 1:L
    parameterValues = get_all_parameter_values(calibmethod, l);

    % Strore parameter Values in structure PV for every time the function is called.
    num = l;
    num = num2str(num);
    eval([ 'PV.parameterValues' num ' =  parameterValues;']);
end

clear l ns parameterValues
delete Store_temp.mat


%% GENERATE MODEL ELEMENTS
% Calls function that calls Dynare to simulate time series and generate 
% model structureswith the parameter values assigned above
ModelStructureM = struct;
ModelStructureoo = struct;
for l = 1:L
    num = num2str(l); % loop identifyer
    save Store_temp l PV dynareFileName num ModelStructureM ModelStructureoo
    % Generate time series
    generate_model_structures(PV, dynareFileName, num);
        
    % Store model structures
    load Store_temp
    eval([ 'ModelStructureM.M_' num ' = M_;']);
    eval([ 'ModelStructureoo.oo_' num ' = oo_;']);
    delete Store_temp
end
tic

%% PERFORM SPECTRAL DECOMPOSITION
% Model structures are used to perform spectral decomposition
load Store_permanent
DecompOutput = [];
for l = 1:L
    num = num2str(l); % loop identifyer
    eval([ 'StructureM_ =  ModelStructureM.M_' num ';']);
    eval([ 'Structureoo_ =  ModelStructureoo.oo_' num ';']);
    % Call function that performs spectral deccomposition
    [decompOutput, decompVars, decompShocks] = ...
            perform_spectral_decomp(StructureM_, Structureoo_);
    DecompOutput = [DecompOutput decompOutput];
end


%% EVALUATE MOMENTS
% Compute median and percentiles over all elements of DecompOutput for all
% loops
Ns = size(decompOutput,2);
MED = []; PER = [];
for j = 1:Ns
     med = median(DecompOutput(:,j:Ns:Ns*L),2);
     per = prctile(DecompOutput(:,j:Ns:Ns*L),[lp,up],2);
     MED = [MED med];   % Includes medians of spectral decomposition for 
                        % all variables (rows) and shocks (columns)
     PER = [PER per];   % Includes percentiles of spectral decomposition 
                        % for all variables (rows) and shocks (columns)
end

disp(decompVars);
disp(decompShocks);
disp(MED);

toc

beep



